#define FANOUT 4

typedef int size_hl;

int child(int p) { return FANOUT * p + 1; };
int parent(int c) { return (c - 1) / FANOUT; };

template<class T, class Pred>
void siftdown(T head, T last, Pred comp, int p) {
  typename Pred::first_argument_type v = head[p];
  typename Pred::first_argument_type max;
  int                                c = child(p);
  int                                maxp;
  const int                          n = last - head;
  const int                          nf = n - FANOUT;

  for ( ; c <= nf; p = c, c = child(p)) {
    maxp = c; max = head[c];
    if (comp(max, head[c + 1])) { maxp = c + 1; max = head[maxp]; };
    if (comp(max, head[c + 2])) { maxp = c + 2; max = head[maxp]; };
#if FANOUT > 3
    if (comp( max, head[c + 3])) { maxp = c + 3; max = head[maxp]; };
#endif
#if FANOUT >= 8
    if (comp(max, head[c + 4])) { maxp = c + 4; max = head[maxp]; };
    if (comp(max, head[c + 5])) { maxp = c + 5; max = head[maxp]; };
    if (comp(max, head[c + 6])) { maxp = c + 6; max = head[maxp]; };
    if (comp(max, head[c + 7])) { maxp = c + 7; max = head[maxp]; };
#endif
#if FANOUT == 16
    if (comp(max, head[c +  8])) { maxp = c +  8; max = head[maxp]; };
    if (comp(max, head[c +  9])) { maxp = c +  9; max = head[maxp]; };
    if (comp(max, head[c + 10])) { maxp = c + 10; max = head[maxp]; };
    if (comp(max, head[c + 11])) { maxp = c + 11; max = head[maxp]; };
    if (comp(max, head[c + 12])) { maxp = c + 12; max = head[maxp]; };
    if (comp(max, head[c + 13])) { maxp = c + 13; max = head[maxp]; };
    if (comp(max, head[c + 14])) { maxp = c + 14; max = head[maxp]; };
    if (comp(max, head[c + 15])) { maxp = c + 15; max = head[maxp]; };
#endif
    c = maxp;
    if (comp(v, max))
	  head[p] = max;
    else {
	  head[p] = v; 
	  return;
	}
  }

  if ( c < n ) {
    maxp = c; max = head[c];
    switch (n - c) {
#if FANOUT == 16
      case 15: if (comp(max, head[c + 14])) { maxp = c + 14; max = head[maxp]; };
      case 14: if (comp(max, head[c + 13])) { maxp = c + 13; max = head[maxp]; };
      case 13: if (comp(max, head[c + 12])) { maxp = c + 12; max = head[maxp]; };
      case 12: if (comp(max, head[c + 11])) { maxp = c + 11; max = head[maxp]; };
      case 11: if (comp(max, head[c + 10])) { maxp = c + 10; max = head[maxp]; };
      case 10: if (comp(max, head[c +  9])) { maxp = c +  9; max = head[maxp]; };
      case  9: if (comp(max, head[c +  8])) { maxp = c +  8; max = head[maxp]; };
      case  8: if (comp(max, head[c +  7])) { maxp = c +  7; max = head[maxp]; };
#endif
#if FANOUT >= 8
	  case 7: if (comp(max, head[c + 6])) { maxp = c + 6; max = head[maxp]; };
      case 6: if (comp(max, head[c + 5])) { maxp = c + 5; max = head[maxp]; };
      case 5: if (comp(max, head[c + 4])) { maxp = c + 4; max = head[maxp]; };
      case 4: if (comp(max, head[c + 3])) { maxp = c + 3; max = head[maxp]; };
#endif
#if FANOUT > 3 
      case 3: if (comp(max, head[c+2])) { maxp = c + 2; max = head[maxp]; };
#endif
      case 2: if (comp(max, head[c+1])) { maxp = c + 1; max = head[maxp]; };
    }
    if (comp(v, max )) {
	  head[p] = max;
	  head[maxp] = v;
	} else
	  head[p] = v;
  } else
	head[p] = v;
};

template <class T, class Pred>
void siftup(T head, T last, Pred comp) {
  typename Pred::first_argument_type t = *(last - 1);
  typename Pred::first_argument_type tt;
  int                                c;
  const int                          n_1less = last - head - 1;

  if (comp(head[0], t)) {
    for (c = n_1less ; c > 0; c = parent(c))
      head[c] = head[parent(c)];
    head[0] = t;
  } else {
    for (c = n_1less; comp(tt = head[parent(c)], t); c = parent(c))
      head[c] = tt;
    head[c] = t;
  }
}

template<class RAI, class Pred> inline
void make_dheap(RAI first, RAI last, Pred Pr) {
  RAI cur;
  
  cur = first + 1;
  while (cur <= last)
    siftup(first, cur++, Pr);
}

template<class RAI, class Pred> inline
void push_dheap(RAI first, RAI last, Pred Pr) {
  siftup(first, last, Pr);
}

template<class RAI, class Pred> inline
void pop_dheap(RAI first, RAI last, Pred Pr) {
  first[0] = *(last - 1);
  siftdown(first, (last - 1), Pr, 0);
}

template <class RAI, class Pred>
void sort_dheap(RAI first, RAI last, Pred Pr) {
  typename Pred::first_argument_type save;
  
  while (last - first > 1) {
    save = *first;
    pop_dheap(first, last, Pr);
    *(--last) = save;
  }
}
